.ipynb .pdf

repository open issue

Binder

Contents

Ozone Monitoring Tutorial

This tutorial provides guided examples of ozone monitoring using data from the Copernicus Atmosphere Monitoring Service (CAMS). It is divided into three parts:

1. View animation of Antarctic ozone hole

2. Calculate the size of the Antarctic ozone hole

3. View vertical profiles of global ozone

It uses CAMS global reanalysis (EAC4) data freely available from the Atmosphere Data Store (ADS)

Run the tutorial via free cloud platforms:

Install and import packages

!pip install cdsapi

Requirement already satisfied: cdsapi in c:\users\cxcs\anaconda3\lib\site-packages (0.5.1)
Requirement already satisfied: tqdm in c:\users\cxcs\anaconda3\lib\site-packages (from cdsapi) (4.62.3)
Requirement already satisfied: requests>=2.5.0 in c:\users\cxcs\anaconda3\lib\site-packages (from cdsapi) (2.26.0)
Requirement already satisfied: charset-normalizer~=2.0.0 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (2.0.4)
Requirement already satisfied: idna<4,>=2.5 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (3.2)
Requirement already satisfied: certifi>=2017.4.17 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (2021.10.8)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (1.26.7)
Requirement already satisfied: colorama in c:\users\cxcs\anaconda3\lib\site-packages (from tqdm->cdsapi) (0.4.4)

# CDS API
import cdsapi

# Libraries for reading and working with multidimensional arrays
import numpy as np
import xarray as xr

# Libraries for plotting and animating data
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import colors
import matplotlib.animation as animation
import cartopy.crs as ccrs
from IPython.display import HTML

# Disable warnings for data download via API
import urllib3 
urllib3.disable_warnings()

from collections import OrderedDict

Download and read and pre-process ozone data

Download data

Copy your API key into the code cell below, replacing ####### with your key. (Remember, to access data from the ADS, you will need first to register/login https://ads.atmosphere.copernicus.eu and obtain an API key from https://ads.atmosphere.copernicus.eu/api-how-to.)

URL = 'https://ads.atmosphere.copernicus.eu/api/v2'

# Replace the hashtags with your key:
KEY = '##########################################'

Here we specify a data directory into which we will download our data and all output files that we will generate:

DATADIR = './'

For this tutorial, we will use CAMS Global Reanalysis (EAC4) data. The code below shows the subset characteristics that we will extract from this dataset as an API request.

Note

Before running this code, ensure that you have accepted the terms and conditions. This is something you only need to do once for each CAMS dataset. You will find the option to do this by selecting the dataset in the ADS, then scrolling to the end of the Download data tab.

c = cdsapi.Client(url=URL, key=KEY)
c.retrieve(
    'cams-global-reanalysis-eac4',
    {
        'variable': 'total_column_ozone',
        'date': '2020-07-01/2021-01-31',
        'time': '00:00',
        'format': 'netcdf',
        'area': [
            0, -180, -90,
            180,
        ],
    },
    f'{DATADIR}/TCO3_202007-202101_SHem.nc')

2022-09-04 23:24:35,351 INFO Welcome to the CDS
2022-09-04 23:24:35,358 INFO Sending request to https://ads.atmosphere.copernicus.eu/api/v2/resources/cams-global-reanalysis-eac4
2022-09-04 23:24:35,438 INFO Request is completed
2022-09-04 23:24:35,444 INFO Downloading https://download-0001-ads-clone.copernicus-climate.eu/cache-compute-0001/cache/data9/adaptor.mars.internal-1661759327.7060935-9333-16-bc780d63-b991-4a5d-8453-84291afb7547.nc to .//TCO3_202007-202101_SHem.nc (23.8M)
2022-09-04 23:24:37,964 INFO Download rate 9.5M/s                                                                      

Result(content_length=24978808,content_type=application/x-netcdf,location=https://download-0001-ads-clone.copernicus-climate.eu/cache-compute-0001/cache/data9/adaptor.mars.internal-1661759327.7060935-9333-16-bc780d63-b991-4a5d-8453-84291afb7547.nc)

Import and inspect data

Here we read the data we downloaded into an xarray Dataset, and view its contents:

fn = f'{DATADIR}/TCO3_202007-202101_SHem.nc'
ds = xr.open_dataset(fn)
ds

<xarray.Dataset>
Dimensions:    (longitude: 480, latitude: 121, time: 215)
Coordinates:
  * longitude  (longitude) float32 -180.0 -179.2 -178.5 ... 177.8 178.5 179.2
  * latitude   (latitude) float32 0.0 -0.75 -1.5 -2.25 ... -88.5 -89.25 -90.0
  * time       (time) datetime64[ns] 2020-07-01 2020-07-02 ... 2021-01-31
Data variables:
    gtco3      (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2022-08-29 07:48:48 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • longitude: 480
  • latitude: 121
  • time: 215
• Coordinates: (3)
  • longitude
    (longitude)
    float32
    -180.0 -179.2 ... 178.5 179.2

    units :

    degrees_east

    long_name :

    longitude

    array([-180.  , -179.25, -178.5 , ...,  177.75,  178.5 ,  179.25],
          dtype=float32)

  • latitude
    (latitude)
    float32
    0.0 -0.75 -1.5 ... -89.25 -90.0

    units :

    degrees_north

    long_name :

    latitude

    array([  0.  ,  -0.75,  -1.5 ,  -2.25,  -3.  ,  -3.75,  -4.5 ,  -5.25,  -6.  ,
            -6.75,  -7.5 ,  -8.25,  -9.  ,  -9.75, -10.5 , -11.25, -12.  , -12.75,
           -13.5 , -14.25, -15.  , -15.75, -16.5 , -17.25, -18.  , -18.75, -19.5 ,
           -20.25, -21.  , -21.75, -22.5 , -23.25, -24.  , -24.75, -25.5 , -26.25,
           -27.  , -27.75, -28.5 , -29.25, -30.  , -30.75, -31.5 , -32.25, -33.  ,
           -33.75, -34.5 , -35.25, -36.  , -36.75, -37.5 , -38.25, -39.  , -39.75,
           -40.5 , -41.25, -42.  , -42.75, -43.5 , -44.25, -45.  , -45.75, -46.5 ,
           -47.25, -48.  , -48.75, -49.5 , -50.25, -51.  , -51.75, -52.5 , -53.25,
           -54.  , -54.75, -55.5 , -56.25, -57.  , -57.75, -58.5 , -59.25, -60.  ,
           -60.75, -61.5 , -62.25, -63.  , -63.75, -64.5 , -65.25, -66.  , -66.75,
           -67.5 , -68.25, -69.  , -69.75, -70.5 , -71.25, -72.  , -72.75, -73.5 ,
           -74.25, -75.  , -75.75, -76.5 , -77.25, -78.  , -78.75, -79.5 , -80.25,
           -81.  , -81.75, -82.5 , -83.25, -84.  , -84.75, -85.5 , -86.25, -87.  ,
           -87.75, -88.5 , -89.25, -90.  ], dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-07-01 ... 2021-01-31

    long_name :

    time

    array(['2020-07-01T00:00:00.000000000', '2020-07-02T00:00:00.000000000',
           '2020-07-03T00:00:00.000000000', ..., '2021-01-29T00:00:00.000000000',
           '2021-01-30T00:00:00.000000000', '2021-01-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • gtco3
    (time, latitude, longitude)
    float32
    ...

    units :

    kg m**-2

    long_name :

    GEMS Total column ozone

    standard_name :

    atmosphere_mass_content_of_ozone

    [12487200 values with dtype=float32]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2022-08-29 07:48:48 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data9/adaptor.mars.internal-1661759327.7060935-9333-16-bc780d63-b991-4a5d-8453-84291afb7547.nc /cache/tmp/bc780d63-b991-4a5d-8453-84291afb7547-adaptor.mars.internal-1661759315.736128-9333-17-tmp.grib

We can see that the data has three coordinate dimensions (longitude, latitude and time) and one variable, total column ozone. By inspecting the coordinates we can see the data is of the southern hemisphere from 1st July 2020 to 31 January 2021 at 00:00 UTC each day. This includes the entire period in which the Antarctic ozone hole appears.

Note

This is only a subset of the available data from the ADS, which includes global data at 3 hourly resolution (in addition to monthly averages), from 2003 to the present, and at 60 model levels (vertical layers) in the atmosphere.

To facilitate further processing, we convert the Xarray Dataset into an Xarray Data Array containing the single variable of total column ozone.

tco3 = ds['gtco3']

Unit conversion

We can see from the attributes of our data are represented as mass concentration of ozone, in units of kg m**-2. We would like to convert this into Dobson Units, which is a standard for ozone measurements. The Dobson unit is defined as the thickness (in units of 10 μm) of a layer of pure gas (in this case O3) which would be formed by the total column amount at standard conditions for temperature and pressure.

convertToDU = 1 / 2.1415e-5

tco3 = tco3 * convertToDU

View southern hemisphere ozone hole

Let us now visualise our data in form of maps and animations.

We will first define the colour scale we would like to represent our data.

Define colour scale

Extract range of values

tco3max = tco3.max()
tco3min = tco3.min()
tco3range = (tco3max - tco3min)/20.

Define colourmap

cmap = plt.cm.jet
norm = colors.BoundaryNorm(np.arange(tco3min, tco3max, tco3range), cmap.N)

Plot map

We will now plot our data. Let us begin with a map for the first time step, 1 July 2020.

fig = plt.figure(figsize=(5, 5)) 
ax = plt.subplot(1,1,1, projection=ccrs.Orthographic(central_latitude=-90)) 
ax.coastlines(color='black') # Add coastlines
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') 
ax.set_title(f'Total Column Ozone, {str(tco3.time[0].values)[:-19]}', fontsize=12) 
im = plt.pcolormesh(tco3.longitude.values, tco3.latitude.values, tco3[0,:,:],
                cmap=cmap, norm=norm, transform=ccrs.PlateCarree())
cbar = plt.colorbar(im,fraction=0.046, pad=0.04) 
cbar.set_label('Total Column Ozone (DU)') 

Create animation

Let us now view the entire time series by means of an animation that shows all time steps. We see from the dataset description above, that the time dimension in our data includes 215 entries. These correspond to 215 time steps, i.e. 215 days, between 1 July 2020 and 31 January 2021. The number of frames of our animation is therefore 215.

frames = 215

Now we will create an animation. The last cell, which creates the animation in Javascript, may take a few minutes to run.

def animate(i):
    array = tco3[i,:,:].values
    im.set_array(array.flatten())
    ax.set_title(f'Total Column Ozone, {str(tco3.time[i].values)[:-19]}', fontsize=12)

ani = animation.FuncAnimation(fig, animate, frames, interval=100)
ani.save(f'{DATADIR}/TCO3_202007-202101_SHem.gif')
HTML(ani.to_jshtml())

2022-09-04 23:27:35,654 INFO Animation.save using <class 'matplotlib.animation.FFMpegWriter'>
2022-09-04 23:27:35,655 INFO MovieWriter._run: running command: ffmpeg -f rawvideo -vcodec rawvideo -s 360x360 -pix_fmt rgba -r 10.0 -loglevel error -i pipe: -filter_complex 'split [a][b];[a] palettegen [p];[b][p] paletteuse' -y .//TCO3_202007-202101_SHem.gif
2022-09-04 23:28:36,187 INFO Animation.save using <class 'matplotlib.animation.HTMLWriter'>

Once Loop Reflect

ACTIVITY:
In what date range does the ozone hole appear greatest? In the next section we will quantify the size of the ozone hole at each time step.

Calculate size of ozone hole

In this section we will calculate the actual size of the ozone hole, and plot its evolution over our time series. We will then apply the exact same method on monthly averaged data over the past decades to compare ozone hole extents at each yearly occurence since 2003.

Define constants

In order to calculate the size of the ozone hole, we need to define some constants. These include the following:

1. Ozone hole limit, i.e. minimum ozone, in Dobson units, beneath which we consider there to be a hole in the ozone layer.

OZONE_HOLE_LIMIT = 220.

1. Radius of the Earth, in order to calculate the size of each geographic grid cell of our data.

Rearth = 6371009.

Calculate area of each geographic grid cell of our data

Here we calculate the size of each grid cell of our data, which varies as a function of latitude. The formula we apply is the following:

area = |(sin(top latitude) - sin(bottom latitude))| * |Δlongitude| * Earth radius squared

First we define a function to apply this formula, taking as arguments the min and max of the latitudes and longitudes represented by each grid cell:

def geo_area_array(lats1,lons1,lats2,lons2):
    area= np.abs(np.sin(np.deg2rad(lats1)) - np.sin(np.deg2rad(lats2))) * np.abs(np.deg2rad(lons1-lons2))
    area = area * Rearth * Rearth 
    return area

Then we apply this function to our data, calculating first the lat/lon extents of each cell (0.75 is the spatial resolution, in degrees):

tco3['areas'] = geo_area_array(tco3['latitude']-0.75/2, tco3['longitude']-0.75/2, 
                             tco3['latitude']+0.75/2, tco3['longitude']+0.75/2)

Mask data belonging to Antarctic ozone hole

We now calculate the grid cells that meet the criteria of belonging to the Antarctic ozone hole: i.e. they have less than the minimum threshold of ozone, and they are below -60 degrees latitude.

Here we define a mask with the threshold conditions for ozone and latitude:

mask = tco3.where((tco3 < OZONE_HOLE_LIMIT) & (tco3["latitude"] < -60.))

All data points that do not meet the conditions are set to NaN (Not a Number) in the resulting array. Dividing this array by itself gives us a mask with each valid data point set to 1. Multiplying this by the corresponding areas leaves us with an array of grid cell areas that meet the conditions of belonging to the Antarctic ozone hole.

area = (mask / mask) * mask['areas']

Calculate ozone hole area

If we now sum these grid cell areas, we have the total ozone hole extent. This is applied across each time step. We multiply by the number 1e-12 to convert from meters squared to million kilometers squared.

ozone_hole = area.sum(dim=["latitude", "longitude"], skipna=True) * 1e-12

Let us update the Data Array attributes:

ozone_hole.attrs['long_name'] = 'Ozone hole area'
ozone_hole.attrs['units'] = 'million km^2'

Plot ozone hole area at each time step

Finally we can plot the evolution of the ozone hole, in million km squared, throughout the time series.

ozone_hole.plot()

[<matplotlib.lines.Line2D at 0x27923c20250>]

We can also calculate the maximum extent of the ozone hole:

index = ozone_hole.argmax() # Here we find the index of the maximum value
print('Max extent of O3 hole:', ozone_hole[index].values, 'million km2')

Max extent of O3 hole: 23.757296 million km2

… and find out when this was reached:

print('This was reached on:', ozone_hole[index].time.values)

This was reached on: 2020-09-28T00:00:00.000000000

ACTIVITY:
Now repeat the steps above with the data below (see API request in cell below). See what you find out about the evolution of the ozone hole since 2003!

c.retrieve(
    'cams-global-reanalysis-eac4-monthly',
    {
        'format': 'netcdf',
        'variable': 'total_column_ozone',
        'year': [
            '2003', '2004', '2005',
            '2006', '2007', '2008',
            '2009', '2010', '2011',
            '2012', '2013', '2014',
            '2015', '2016', '2017',
            '2018', '2019', '2020',
            '2021',
        ],
        'month': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
        ],
        'product_type': 'monthly_mean',
    },
    f'{DATADIR}/TCO3_monthly_2003-2021.nc')

Ozone vertical profile

So far we have concentrated on total column ozone. In this final part of the tutorial we will look at vertical profiles of ozone to see how ozone concentration differs across various altitudes.

This requires a new data request, including multi level data at 25 pressure levels represented in hectopascals (hPa).

Note

The ADS provides freely available data at 60 atmospheric levels (model levels), but most of this is on slower access tapes. For the sake of simplicity, we will use in this tutorial only data from the 25 pressure levels.

c.retrieve(
    'cams-global-reanalysis-eac4-monthly',
    {
        'variable': 'ozone',
        'pressure_level': [
            '1', '2', '3',
            '5', '7', '10',
            '20', '30', '50',
            '70', '100', '150',
            '200', '250', '300',
            '400', '500', '600',
            '700', '800', '850',
            '900', '925', '950',
            '1000',
        ],
        'year': '2021',
        'month': '07',
        'product_type': 'monthly_mean',
        'format': 'netcdf',
    },
    f'{DATADIR}/O3_2021-07.nc')

2022-08-31 22:59:31,303 INFO Welcome to the CDS
2022-08-31 22:59:31,309 INFO Sending request to https://ads.atmosphere.copernicus.eu/api/v2/resources/cams-global-reanalysis-eac4-monthly
2022-08-31 22:59:31,387 INFO Request is queued
2022-08-31 22:59:32,464 INFO Request is running
2022-08-31 22:59:34,046 INFO Request is completed
2022-08-31 22:59:34,060 INFO Downloading https://download-0001-ads-clone.copernicus-climate.eu/cache-compute-0001/cache/data6/adaptor.mars.internal-1661983111.4206462-29418-12-b8063567-b4a8-48a6-873e-e21165220e7d.nc to ./DATA/O3_202107.nc (5.5M)
2022-08-31 22:59:35,065 INFO Download rate 5.5M/s                                                                      

Result(content_length=5788292,content_type=application/x-netcdf,location=https://download-0001-ads-clone.copernicus-climate.eu/cache-compute-0001/cache/data6/adaptor.mars.internal-1661983111.4206462-29418-12-b8063567-b4a8-48a6-873e-e21165220e7d.nc)

fn = f'{DATADIR}/O3_2021-07.nc'
ds = xr.open_dataset(fn)
ds

<xarray.Dataset>
Dimensions:    (longitude: 480, latitude: 241, level: 25, time: 1)
Coordinates:
  * longitude  (longitude) float32 0.0 0.75 1.5 2.25 ... 357.0 357.8 358.5 359.2
  * latitude   (latitude) float32 90.0 89.25 88.5 87.75 ... -88.5 -89.25 -90.0
  * level      (level) int32 1 2 3 5 7 10 20 30 ... 700 800 850 900 925 950 1000
  * time       (time) datetime64[ns] 2021-07-01
Data variables:
    go3        (time, level, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2022-08-31 21:58:31 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

xarray.Dataset

• Dimensions:
  • longitude: 480
  • latitude: 241
  • level: 25
  • time: 1
• Coordinates: (4)
  • longitude
    (longitude)
    float32
    0.0 0.75 1.5 ... 357.8 358.5 359.2

    units :

    degrees_east

    long_name :

    longitude

    array([  0.  ,   0.75,   1.5 , ..., 357.75, 358.5 , 359.25], dtype=float32)

  • latitude
    (latitude)
    float32
    90.0 89.25 88.5 ... -89.25 -90.0

    units :

    degrees_north

    long_name :

    latitude

    array([ 90.  ,  89.25,  88.5 , ..., -88.5 , -89.25, -90.  ], dtype=float32)

  • level
    (level)
    int32
    1 2 3 5 7 ... 850 900 925 950 1000

    units :

    millibars

    long_name :

    pressure_level

    array([   1,    2,    3,    5,    7,   10,   20,   30,   50,   70,  100,  150,
            200,  250,  300,  400,  500,  600,  700,  800,  850,  900,  925,  950,
           1000])

  • time
    (time)
    datetime64[ns]
    2021-07-01

    long_name :

    time

    array(['2021-07-01T00:00:00.000000000'], dtype='datetime64[ns]')

• Data variables: (1)
  • go3
    (time, level, latitude, longitude)
    float32
    ...

    units :

    kg kg**-1

    long_name :

    Ozone mass mixing ratio (full chemistry scheme)

    standard_name :

    mass_fraction_of_ozone_in_air

    [2892000 values with dtype=float32]

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2022-08-31 21:58:31 GMT by grib_to_netcdf-2.25.1: /opt/ecmwf/mars-client/bin/grib_to_netcdf.bin -S param -o /cache/data6/adaptor.mars.internal-1661983111.4206462-29418-12-b8063567-b4a8-48a6-873e-e21165220e7d.nc /cache/tmp/b8063567-b4a8-48a6-873e-e21165220e7d-adaptor.mars.internal-1661983109.786136-29418-20-tmp.grib

Pressure levels in logarithmic scale

To facilitate visualisation of this data, we will create a new coordinate of the array with the pressure levels converted into logarithms.

ds['level_log10'] = np.log10(ds['level'])

ds = ds.set_coords('level_log10')

Reorganise data and convert units

Here we create an Xarray Data Array object from the Dataset and remove the dimension of time, which has only one entry.

da = ds['go3']

da = da.squeeze('time')

We will now convert the units of our data from mass mixing ratio (kg/kg) to Volume Mixing Ratio (VMR), in units of parts per billion (ppbv). We will also update the attributes of our array to reflect this.

da = 28.9644 / 47.9982 * 1e9 * da

da.attrs['long_name'] = 'Ozone Volume Mixing Ratio'
da.attrs['units'] = 'ppbv'

Convert longitudes to +/- 180 degrees

Below, we would like to plot a 2 dimensional map of our data. Note however that the longitudinal coordinates are given in the range 0 to 360 degrees. To facilitate data visualisation we will convert the longitudes into a +/- 180 degree grid.

da = da.assign_coords(longitude=(((da.longitude + 180) % 360) - 180)).sortby('longitude')

Visualise geographical distribution of ozone at each layer

We will now visualise the global distribution of ozone at particular layers of the atmosphere. The description of our data above shows us that we have 25 layers of the atmosphere represented in our data, from 1000 hPa (lowest layer of atmosphere) to 1 (uppermost layer). These can be indexed in the first dimension of our array, where the 24th index corresponds to the lowest layer of the atmosphere, and the 0th represents the highest.

da[24,:,:].plot()

<matplotlib.collections.QuadMesh at 0x2bfb3878100>

ACTIVITY:
View the ozone at different levels of the atmosphere. Note the difference in the scale on the right!

Spatial aggregation

We will now visualise the global average of ozone at each layer in a profile plot. To do this, we must first average over the longitudinal and latitudinal dimensions to obtain a single figure at each layer. When averaging over the latitudinal dimension we need to account for the variation in area represented by each latitudinal grid cell. We do this by applying weights corresponding to the cosine of the latitudes.

weights = np.cos(np.deg2rad(da.latitude))
weights.name = "weights"
weighted = da.weighted(weights)

o3 = weighted.mean(dim=["latitude", "longitude"])

Plot ozone profile

We can now create our profile plot of the vertical distribution of ozone.

fig, ax = plt.subplots(1, 1, figsize = (6, 9))

ax.set_title('Ozone global average July 2021', fontsize=12)
ax.set_ylabel('Log pressure (10^x hPa)')
ax.set_xlabel('O3 (ppbv)')
ax.invert_yaxis()
ax.plot(o3, o3.level_log10)

fig.savefig(f'{DATADIR}/O3_2021-07_profile.png')

ACTIVITY:
At what level of the atmosphere is the concentration of ozone the greatest? Hint: To find the index of the maximum value of the array, use the same .argmax() function you used when calculating the size of the ozone hole above!

previous

European Air Quality Index Calculation